#ifndef _accumulate_hpp_
#define _accumulate_hpp_

#include <numeric>

#include "vector_traits.hpp"

namespace core {

    namespace math {

        namespace detail {

            template<class V, class F, std::size_t n>
                struct accumulate_unroll_loop
                {
                    typedef typename field<V>::type  field_type;
                    inline static field_type eval(const V& v, const field_type& init_value, F binary_op)
                    {
                        return binary_op(accumulate_unroll_loop<V,F,n-1>::eval(v,init_value,binary_op),v[n]);
                    }
                };

            // break condition
            template<class V, class F>
                struct accumulate_unroll_loop<V,F,0>
                {
                    typedef typename field<V>::type  field_type;
                    inline static field_type eval(const V& v, const field_type& init_value, F binary_op)
                    {
                        return binary_op(init_value,v[0]);
                    }
                };

            struct accumulate_static_impl
            {
                template<class Vector, class BinaryFunctor>
                    inline static typename field<Vector>::type eval(const Vector& v, const typename field<Vector>::type& init_value, BinaryFunctor binary_op)
                    {
                        typedef typename field<Vector>::type  field_type;
                        return accumulate_unroll_loop<Vector,BinaryFunctor,
                               dimension<Vector>::value-1 >::eval(v,init_value,binary_op);
                    }
            };

            struct accumulate_dynamic_impl
            {
                template<class Vector, class BinaryFunctor>
                    inline static typename field<Vector>::type eval(const Vector& v, const typename field<Vector>::type& init_value, BinaryFunctor binary_op)
                    {
                        typedef typename field<Vector>::type  scalar_type;
                        using std::accumulate;
                        return accumulate(v.begin(),v.end(),init_value,binary_op);
                    }
            };

            struct accumulate_scalar_impl
            {
                template<class S, class F>
                    inline static S eval(const S& value, const S& init_value, F binary_op)
                    {
                        return binary_op(init_value,value);
                    }
            };

            template<class T, class F>
                struct choose_accumulate_impl
                {
                    //check if this is a vector with static dimension
                    typedef typename boost::mpl::eval_if<is_scalar<T>,
                                                boost::mpl::identity<accumulate_scalar_impl>,
                                                boost::mpl::eval_if_c<(dimension<T>::value>0),
                                                                    boost::mpl::identity<accumulate_static_impl>,
                                                                    boost::mpl::identity<accumulate_dynamic_impl>
                                                                    >
                                                >::type  impl_type;

                    typedef typename field<T>::type scalar_type;
                    inline static scalar_type eval(const T& v, const scalar_type& init_value, F binary_op)
                    {
                        return impl_type::template eval<T,F>(v,init_value,binary_op);
                    }
                };

        } // namespace detail

        /** \brief std::accumulate-like method with optional loop unrolling for static arrays and vectors
         *
         */
        template<class Vector, class BinaryFunctor>
            typename field<Vector>::type accumulate(const Vector& v, const typename field<Vector>::type& init_value, BinaryFunctor& binary_op)
            {
                typedef typename detail::choose_accumulate_impl<Vector,BinaryFunctor&> impl_type;
                return impl_type::eval(v,init_value,binary_op);
            }

        template<class Vector, class BinaryFunctor>
            typename field<Vector>::type accumulate(const Vector& v, const typename field<Vector>::type& init_value, const BinaryFunctor& i_binary_op)
            {
                BinaryFunctor binary_op(i_binary_op);
                typedef typename detail::choose_accumulate_impl<Vector,BinaryFunctor&> impl_type;
                return impl_type::eval(v,init_value,binary_op);
            }

    } // namespace math

} // namespace core

#endif // _accumulate_hpp_
